Draft version July 9, 2009 

Preprint typeset using I4TgX style emulateapj v. 03/07/07 



OS 

o 
o 

(N 



oo 

S 1 

6 
& 

CO 

> 
O 
00 

in 



oo 
o 



THE BLAZAR SEQUENCE AND THE COSMIC GAMMA-RAY BACKGROUND RADIATION 

IN THE FERMI ERA 

YOSHIYUKI INOUE AND TOMONORI TOTANI 
Department of Astronomy, Kyoto University, Kitashirakawa, Sakyo-ku, Kyoto 606-8502, Japan 

Draft version July 9, 2009 

ABSTRACT 

We present a new model of the blazar gamma-ray luminosity function (GLF) and the spectrum of the ex- 
tragalactic gamma-ray background (EGRB), which is consistent with the observed distributions of EGRET 
blazars. The unified sequence of blazar spectral energy distribution (SED) is taken into account to make a non- 
trivial prediction for the EGRB spectrum and more realistic comparison with the data than previous studies. 
We then try to explain the EGRB data by the two AGN populations: one is blazars, and the other is non-blazar 
AGNs that are responsible for the EGRB in the MeV band. We find that ~80% of the EGRB photon flux at 
> 100 MeV can be explained by the sum of the two populations, while ^45 % can be accounted for only by 
blazars. The predicted EGRB spectrum is in agreement with a wide range of the observed data from X-ray 
to GeV, within the systematic uncertainties in the EGRB determination by EGRET. These results indicate that 
AGNs including blazars are the primary source of EGRB. Blazars are dominant in EGRB at higher energy 
bands of > 100 MeV, while non-blazar AGNs dominate at < 100 MeV. Almost all of the EGRB flux from 
blazars will be resolved into discrete sources by the Fermi Gamma-ray Space Telescope, while that from non- 
blazar AGNs will largely remain unresolved. Therefore, comparison between the integrated source counts and 
diffuse EGRB flux as a function of photon energy will give a simple and clear test of our model. Various 
quantitative predictions for Fermi observations are also made. Especially, our model predicts 600-1200 blazars 
in all sky down to 2 x 10~ 9 photons cm" 2 s _1 (> 100 MeV), which is considerably smaller than most of previous 
studies. We find that the fraction of EGRB energy flux absorbed in intergalactic medium (IGM) is not large, 
and the cascade component reprocessed in IGM does not significantly alter the EGRB spectrum. 
Subject headings: cosmology: diffuse radiation - galaxies : active - gamma rays : theory 



1. INTRODUCTION 

The origin of the extragalactic diffuse gamma-ray back- 
ground (EGRB) has been discussed in astrophysics for a long 
time. EGRB was first discovered by the SAS 2 satellite (Fich- 
tel, Simpson, & Thompson 1978; Thompson & Fichtel 1982). 
Subsequently, the EGRB spectrum was confirmed up to ~ 
50 GeV by EGRET (Energetic Gamma-Ray Experiment Tele- 
scope) on board the Compton Gamma Ray Observatory. The 
EGRB flux is about lx 10~ 5 photons cm -2 s" 1 sr" 1 above 
100 MeV and the spectrum is approximately a power law 
with a photon index of ~ -2 in a wide range of ~30 MeV 
- 100 GeV (Sreekumar et al. 1998; Strong, Moskalenko, 
& Reimer 2004a). It should be noted that, however, mea- 
surement of EGRB is not an easy task. The Galactic dif- 
fuse background from cosmic -ray interaction in the Galac- 
tic disk is a strong foreground emission and must be sub- 
tracted to estimate EGRB. Therefore modeling of this fore- 
ground component could induce significant systematic uncer- 
tainties in EGRB measurements (Keshet, Waxman, & Loeb 
2004; Strong, Moskalenko, & Reimer 2004a,b; Kamae, Abe, 
& Koi 2005; Kamae et al. 2006). A possible systematic er- 
ror in the calibration of the EGRET detector may also have 
affected the EGRB determination (Stecker et al. 2008). 

Although several sources of gamma-rays (e.g., clusters of 
galaxies or dark matter annihilation) have been proposed as a 
significant contributor to EGRB [see, e.g., Narumoto & Totani 
(2006) and references therein], active galactic nuclei (AGNs) 
of the blazar class have been thought as the primary candidate 
for the origin of EGRB, since almost all of the extragalac- 
tic gamma-ray sources detected by EGRET are blazars. The 
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blazar contribution to EGRB has been estimated by a number 
of papers (Padovani et al. 1993; Stecker, Salamon, & Malkan 
1993; Salamon & Stecker 1994; Chiang et al. 1995; Stecker 
& Salamon 1996; Chiang & Mukherjee 1998; Mukherjee 
& Chiang 1999; Miicke & Pohl 2000; Narumoto & Totani 
2006; Giommi et al. 2006; Dermer 2007; Pavlidou & Ven- 
ters 2008; Kneiske & Mannheim 2008; Bhattacharya et al. 
2008). These previous studies have derived different results 
by different approaches about the contribution of blazars to 
EGRB, ranging from 20% to 100 %. A brief summary of 
these studi es an d comparison with our new result will be pre- 
sented in §17.11 In our previous study (Narumoto & Totani 
2006, hereafter NT06), we constructed a gamma-ray luminos- 
ity function (GLF) model based on the picture of luminosity- 
dependent density evolution (LDDE), which has been known 
to describe well the evolution of X-ray luminosity function 
(XLF) of AGNs. NT06 found that the LDDE model fits to the 
EGRET blazar data better than the pure-luminosity-evolution 
(PLE) models employed in most of the previous studies. It is 
then found that only 25-50% of EGRB can be explained by 
blazars, with the GLF model parameters that are consistent 
with the EGRET blazar data (flux and redshift distributions). 

In most of past studies including NT06, however, blazar 
spectral energy distributions (SEDs) were assumed to be a 
single or broken power-law for all blazars. In such a mod- 
eling, the predicted EGRB spectrum is almost obviously de- 
termined by the assumed power-law indices. Blazar SED 
has been theoretically and observationally studied in detail 
(Ghisellini & Tavecchio 2008b, and references therein). From 
the theoretical point of view, blazar emission is widely be- 
lieved to be the sum of the synchrotron (radio to UV bands) 
and inverse Compton (dominant in gamma-ray bands) compo- 
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nents produced by the same nonthermal electron population 
accelerated in relativistic jets. The source of target photons 
for the inverse-Compton (IC) component could be either the 
synchrotron photons produced in the jet itself (synchrotron- 
self-Compton, SSC), or photons emitted from the accretion 
disk (external Compton, EC). Multi-wavelength observational 
studies from radio to 7-ray bands have indicated an interest- 
ing feature in blazar SEDs; the synchrotron and Compton 
peak photon energies decrease as the bolometric luminosity 
increases (Fossati et al. 1997, 1998; Ghisellini et al. 1998, 
Donato et al. 2001, Ghisellini & Tavecchio 2008b, Maraschi 
et al. 2008). This is often called as the blazar SED sequence. 
Although the validity of the blazar sequence is currently still 
a matter of debate (e.g., Padovani et al. 2007; Ghisellini & 
Tavecchio 2008b; Maraschi et al. 2008), one can make a 
non-trivial prediction of the EGRB spectrum if this blazar se- 
quence is assumed. In other words, the blazar SED sequence 
can be tested by comparing with the observed EGRB spec- 
trum. 

In this paper we calculate the EGRB flux and spectrum from 
blazars, by constructing a blazar GLF model that is consistent 
with the flux and redshift distributions of the EGRET blazars, 
based on the LDDE scheme and the blazar SED sequence. By 
introducing the blazar SED sequence, we can make a reason- 
able and non-trivial prediction of the EGRB spectrum for the 
first time, which can be compared with the observed EGRB 
spectrum. 

Recently, Inoue, Totani, & Ueda (2008) has showed that 
EGRB in the MeV band can naturally be explained by normal 
(i.e., non-blazar) AGNs that compose the cosmic X-ray back- 
ground. This MeV background component extends to ~ 100 
MeV with a photon index of about 2.8, by the Comptoniza- 
tion photons produced by nonthermal electrons in hot coro- 
nae. Therefore, it should also contribute to the EGRB at < 1 
GeV We will investigate how much fraction of the observed 
EGRB can be explained by the sum of the two components, 
i.e., non-blazar AGNs (dominant at <100 MeV) and blazars 
(dominant at > 100 MeV). 

We will then make quantitative predictions for the Fermi 
Gamma-ray Space Telescope (Atwood et al. 2009; formerly 
known as GLAST) that has successfully been launched on 
June 1 1th, 2008, so that our model can be tested quantitatively 
by the observational data coming in the very near future 1 . 

This paper is organized as follows. In ^2 we will present 
the treatment of the blazar SED sequence in our calculation. 
We will describe the formulations for the blazar GLF model 
taking into account the SED sequence in 5j3j an d me GLF 
parameter determination by fitting to the EGRET data in 50] 
The EGRB will be calculated and compared with the ob- 
served data in $5] and we will present predictions for Fermi 
in ^6] Discussions on a few issues will be given in $2l in- 
cluding a comparison of our results with those in previous 
studies. Conclusions will then be given in ^8] Throughout 
this paper, we adopt the standard cosmological parameters of 
(/7,f! M ,fi A )=(0.7,0.3,0.7). 

2. THE BLAZAR POPULATION AND THE SED SEQUENCE 

1 After we submit the first version of this paper larXiv:0810.3580vl), the 
first Fermi catalog for bright gamma-ray sources including AGNs has appear 
(Abdo et al. 2009a, b). The sample size of blazars is about 100, which is 
bigger than the EGRET catalog by a modest factor of about 2. We confine 
ourselves using only the EGRET data in this work, as the theoretical predic- 
tion in the pie-Fermi era. A much larger number of blazars should be detected 
in the future Fermi sample, which should be compared with our predictions. 



Blazars are defined as the combined class of the two pop- 
ulations of AGNs: BL Lac objects and flat spectrum radio 
quasars (FSRQs) (Urry & Padovani 1995) 2 . BL Lacs are de- 
fined as AGNs whose nuclear non-thermal continuum emis- 
sion is so strong that the rest-frame equivalent width (EW) 
of the strongest optical emission line is narrower than 5 A. 
Such domination of continuum is interpreted as the jet be- 
ing directed towards the observer; the continuum from jet 
is extremely enhanced by the relativistic beaming compared 
with more isotropic line emission from the regions around 
accretion disks. FSRQs are AGNs with emission line EW 
greater than 5 A, but whose spectral index in radio band is 
less than a r < 0.5, where a r is defined by f u cx v~ a ' (see 
Urry & Padovani 1995 for detailed reviews). Although FS- 
RQs show discernible emission lines, the flat radio spectrum 
is also an evidence of the dominating jet component. In fact 
FSRQs show similar properties to BL Lacs, such as rapid 
time variability and strong polarization. These properties in- 
dicate that BL Lacs and FSRQs are physically similar popu- 
lations, and hence they are often combined into a single class 
of blazars. The parameter that discriminates BL Lacs and 
FSRQs is likely to be the luminosity; generally FSRQs have 
larger luminosities than BL Lacs (Fossati et al. 1998, Ghis- 
ellini et al. 1998) 3 . In this paper we also treat the BL Lacs 
and FSRQs as the single population of blazars 4 . 

Fossati et al. (1997, 1998) and Donato et al. (2001, here- 
after D01) constructed an empirical blazar SED model to de- 
scribe the SED sequence, based on fittings to observed SEDs 
from radio to 7-ray bands. These models are comprised of 
the two components (synchrotron and IC), and each of the 
two is described by a linear curve at low photon energies and 
a parabolic curve at high energies, in the plane of log 10 (i/L„) 
and log 1() v. 

Here we construct our own SED sequence model mainly 
based on the SED model of D01, because there is a math- 
ematical discontinuity in the original D01 model. In the 
D01 model, two different mathematical fitting formulae are 
used below and above the luminosity vh v = 10 43 erg s" 1 at 
5 GHz, and the luminosity of the inverse-Compton compo- 
nent suddenly changes with a jump at this critical luminos- 
ity. Our own SED sequence formula is described in Ap- 
pendix in detail, and the discontinuity is avoided there. Once 
the blazar luminosity is specified at a reference frequency 
(e.g., 5 GHz in radio band), this empirical model gives the 
full blazar SED from radio to gamma-ray bands. In FigQ] 
we show this empirical blazar SED sequence model in com- 
parison with the observed SED data (Fossati et al. 1998; 
D01). It should be noted that the observed data are means 
of many blazars in a certain luminosity range, and there may 
be scatter of individual blazar SEDs from the sequence. When 
the bolometric blazar luminosity P = j vh v is specified, the 
blazar luminosity per unit frequency, L v (u,P), is determined 
for any photon frequency v by the blazar sequence model. 
Note that, although blazars must be strongly anisotropic emit- 

2 Sometimes FSRQs are also called as optically violent variable (OVV) 
quasars or highly polarized quasars (HPQs), but they are essentially the same 
populations. 

3 A recently popular classification is that blazars are divided into quasar- 
hosted blazars (QHB) that are the same population as FSRQ, low frequency 
BL Lacs (LBLs), and high frequency BL Lacs (HBLs), in the decreasing 
order of the absolute luminosity (Kubo et al. 1998). 

4 After the submission of the first version of our paper 
larXiv:0810.3580vl), Ghisellini et al. (2009) has shown that this treat- 
ment is adequate by using the latest data of Fermi. 
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FIG. 1 . — The blazar SED sequence. The data points are the average SED of the blazars studied by Fossati et al. (1998) and D01. The solid curves are the 
empirical SED sequence models constructed and used in this paper. The model curves corresponds to the bolometric luminosities of log 10 (P/erg s -1 ) = 49.50, 
48.64, 47.67, 46.37, and 45.99 (from top to bottom). 



ters, L v is defined as the isotropic-equivalent luminosity and 
hence the observed energy flux per unit frequency is given by 
F v [v/(1 +z)] = (1 + z)L„(v 1 P)/[4ird L (z) 2 ], where z is redshift 
and cIl is the standard luminosity distance. 

3. THE MODEL OF GAMMA-RAY LUMINOSITY FUNCTION OF 
BLAZARS 

3.1. Relating Jet Power and X-ray Luminosity ofAGNs 

The cosmological evolution of X-ray luminosity function 
of AGNs has been investigated intensively [Ueda et al. 2003 
(hereafter U03); Hasinger, Miyaji, & Schmidt 2005 (here- 
after H05); Sazonov et al. 2007; Gilli, Comastri, & Hasinger 
2007]. These studies revealed that AGN XLF is well de- 
scribed by the LDDE model, in which peak redshift of density 
evolution increases with AGN luminosity. Here we construct 
two models of blazar GLF based on the two XLFs derived by 
U03 (in hard X-ray band) and H05 (in soft X-ray band), both 
of which are based on the LDDE scheme. The use of LDDE 
in blazar GLF has been supported from the EGRET blazar 
data, because NT06 found that the EGRET data agrees with 
the LDDE model better than PLE models. However, the va- 
lidity of this assumption from a theoretical viewpoint should 
also be examined. 

It is known that radio jet emission is linked to the dissipa- 
tion process in the accretion disk (Falcke & Biermann 1995; 
Falcke et al. 1995; Merloni et al. 2003; Falcke et al. 2004; 
Kording et al.2006). Therefore, it is natural to expect that 
power of blazar jet is correlated with mass accretion rate onto 
super massive black holes, which is also correlated with the 
X-ray luminosity from accretion disk. Therefore we simply 
assume that the bolometric luminosity of radiation from jet, 



P, is proportional to disk X-ray luminosity, Lx- It should be 
noted that Lx is not the observed X-ray luminosity of a blazar 
having a jet luminosity P; when an AGN is observed as a 
blazar (i.e., the jet directed to an observer), its X-ray flux is 
dominated by the radiation from the jet that is much brighter 
than that from the accretion disk. Rather, Lx is the luminosity 
that would be observed from a direction different from the jet. 

A constant P/Lx ratio is realized when, e.g., the jet kinetic 
luminosity (Pk) is efficiently dissipated into blazar bolomet- 
ric luminosity (P), and both P k and L x are proportional to 
the mass accretion rate (m). One should, however, be care- 
ful about the latter condition. Recent observations of X-ray 
binaries indicate that P k is generally proportional to m, but Lx 
is not, when the accretion rate is much lower than the Edding- 
ton limit (i.e., low Eddington ratio) (Gallo et al. 2003; Gallo 
et al. 2005). Such accretion disks are described by radiatively 
inefficient accretion flows (RIAF) rather than the standard ac- 
cretion disk, and Lx is roughly proportional to m 2 in the RIAF 
regime (Kato et al. 1998; Narayan & Quataert 2005). The 
RIAF picture is well consistent also with the accretion flow 
onto the supermassive black hole of the Galaxy (i.e., Sgr A*) 
(see, e.g., Totani 2006 and references therein). 

Accretion rates of X-ray bright AGNs used to derive the 
AGN XLF are generally close to the Eddington limit, oth- 
erwise they are hardly detected by X-ray observations be- 
cause of the rapid decrease of X-ray luminosity with decreas- 
ing Eddington ratio. Therefore, our blazar GLF should be 
considered as that for high-Eddington-ratio AGNs, and low- 
Eddington-ratio AGNs in the RIAF regime could be missed 
in our analysis. Such a low-Eddington-ratio population might 
have a significant contribution to blazar GLF, because we ex- 
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TABLE 1 
The parameters of the AGN XLF 



Ueda et al. 2003 Hasinger et al. 2005 



Ax" 


5.04 X 10~* 


2.62 x 10~ 7 




43 94-H).2i 


43.94±0.11 


71 


0.86±0.15 


0.87 ±0.10 


72 


2.23 ±0.13 


2.57 ±0.16 


z? 


1.9 C 


1.96±0.15 




44.6 C 


44.67 c 


a 


0.335 ±0.07 


0.21 ±0.04 


"> 

Pi 


4.23 ±0.39 


4.7 ±0.3 


-1.5 C 


-1.5±0.7 


ft 


0.0 rf 


0.7 ±0.3 


02 


0.0 rf 


0.6±0.8 


a In units of Mpc 


-3 




In units of ergs 


s-' 1 . 





c These quantities are treated as fixed parameters in each XLF model. 



The indices j3\, $2 are treated as constants in U03. 



pect Pk oc m in contrast to Lx oc m 2 in the RIAF mode. How- 
ever, the black hole mass function predicted by time inte- 
gration of X-ray AGN LF is consistent with the local black 
hole mass function inferred from the black-hole-mass versus 
bulge-mass relation, indicating that black hole mass grows 
mainly in the high Eddington ratio phase (e.g., Marconi et al. 
2004). If this is correct, cosmic background radiation from 
jet activities should also be dominated by AGNs in high Ed- 
dington ratio phase. Hence, it is reasonable to expect that a 
significant part of EGRB flux can be accounted for by blazars 
with high Eddington ratio phase, whose GLF evolution is de- 
scribed by LDDE. 

3.2. Model Formulations 

In this paper we describe the blazar GLF in terms of vL v 
luminosity at a reference rest-frame photon energy e le f, res = 
100 MeV, i.e., L 7 = (e ref:KS /h p ) L„(e letle Jh p ,P), where h p is 
the Planck constant. According to the assumption justified in 
33. 11 we simply relate the bolometric blazar luminosity P and 
disk X-ray luminosity by the parameter q, as: 



P=\Q q L x 



(1) 



Here, we define the disk luminosity Lx to be that in the rest- 
frame 2-10 and 0.5-2 keV bands for the hard XLF (U03) and 
the soft XLF (H05), respectively. Thus, L 7 and Lx have been 
related through P. 
The blazar GLF p 7 is then obtained from the AGN XLF px, 

as 

dL x 

|0 7 (L 7 ,z) = n-^-px(Lx,z), (2) 

where p 1 and px are the comoving number densities per unit 
gamma-ray and X-ray luminosity, respectively. The param- 
eter k is a normalization factor, representing the fraction of 
AGNs observed as blazars. In our GLF model, we adopt the 
same form in U03 and H05 for p x , as: 



Px (Lx , z) = px (Lx , 0)f(L x , z) , 



(3) 



where px(Lx,0) is the AGN XLF at present. This is charac- 
terized by the faint-end slope index 71 , the bright-end slope 
index 72, and the break luminosity L x , as: 



Px(L x ,0): 



L x ln(10) 



(4) 



where Ax is the normalization parameter 5 having a dimension 
of volume" 1 . The function f{Lx,z) describes the density evo- 
lution, which is given by the following forms: 



f(L x ,z) = 



f(LxMLx)) (1 



1+2 



l'2 



z < Zc(Lx), 
' z > Zc(Lx), 



where z c is the redshift of evolutionary peak, given as 



Zc(Lx)-- 



z* L x > L a , 

z* c (L x /L a ) a L x <L a , 



(5) 



(6) 



The evolutionary indices p\ and p2 are described by using the 
parameters p\, p%, Pi, and P2: 

p 1 =p* l +f3 l (log l0 L x -44.0), (7) 
P2=P* 2 +P2(log w L x -44.0). (8) 

The parameters obtained by the fit to the observed data of X- 
ray AGNs in U03 and H05 are shown in Table [Q 

When we calculate the EGRB flux, it diverges if 71 > 1 
and GLF is integrated down to L 1 — > 0. Therefore we set the 
minimum of the gamma-ray luminosity as L 7iimn =10 43 erg s" 1 
in the EGRB calculation, because we will find that there is 
no EGRET blazars below this value (see Fig. However, 
it should be kept in mind that there might be a considerable 
contribution to the EGRB flux from blazars below L 7 . m i n when 
71 > I- 

4. GAMMA-RAY LUMINOSITY FUNCTION DETERMINED BY THE 
EGRET BLAZAR DATA 

4.1. The Maximum Likelihood Method 

We use the maximum likelihood method to search for the 
best-fit model parameters of the blazar GLF to the distri- 
butions of the observed quantities of the EGRET blazars 
(gamma-ray flux and redshift). The analysis method and the 
data used are essentially the same as those in NT06. 

Observed gamma-ray flux F- 1 of EGRET blazars are photon 
flux at e 7 > e m ; n . bs = 100 MeV in photons cm" 2 s" 1 , where 
emin,obs is in the observer's frame. For a given redshift, F 7 can 
be related to P by the blazar SED sequence model as follows: 



l+Z 



4nd L (zy 



n,obs 



h p v 



dv 



(9) 



Since P has been related to the gamma-ray luminosity L 7 by 
the SED sequence, one can calculate L 1 (z,F J ) for a given set 
of z and F y through P. 

A specified GLF model predicts the distribution function 
d 3 N/dzdF 7 dil of redshift z, flux F 1 , and the location of a 
blazar in the sky specified by a solid angle f2, which is given 
as: 



d 3 N(z,F 1 ,Q) dV 



dz 



dzdF 7 dfl 



P 1 (L 1 ,z)e(F 1 ,z) 



(10) 



x9LF 7 -F 7iUm (0)J, 



where dV / dz is the comoving volume element per unit solid 
angle, the step function (0(x) = 1 and for x > and x < 0, 
respectively), and Fj,\i m (il) the sensitivity limit of EGRET for 
point sources that is a function of the Galactic latitude. The 
functional form of f 7 .ii m (fi) is given in NT06. The detection 
efficiency e(F 7 ,z) represents the identification probability as 

5 The factor of (Lx In 10)~' appeal's because Ax is defined as the pre-factor 
of the number density per unit log 10 Lx. 
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TABLE 2 

Best- fit parameters for blazar GLF 





U03(4) 


U03( ? , 7l ) 


H05(?) 


H05( 9 , 7l ) 


71 

K 


4 Q9+0.21 
4 - yz -o.io 
0.86" 

1.7 x icr 6 


4 01+0.25 

0.93 ±0.1 3 
1.5 x 10~ 6 


c 9O+0.26 
J -^'-0.22 

0.87" 
9.5 x 10- 6 


5 35 +<) - 25 
lli -0.12 

6.0 x 10" 6 


KS test probabilities 


z 


53.8% 


86.9% 


0.15% 


33.0% 




25.6% 


48.4% 


0.08% 


28.5% 



NOTE. — The best-fit values of the model parameters (q, 71 , k) obtained from the 
maximum likelihood analysis. The KS probabilities of the best-fit models are also shown 
for the redshift and gamma-ray luminosity distributions in the last two rows. 
a These parameters are fixed at the original AGN XLF values in these analysis, and the 
original values are shown. 



a blazar by finding a radio counterpart. This is described §2.3 
in NT06, but here we modified the relation between the radio 
luminosity L R and L 7 from a simple linear relation in NT06 
to that predicted from our blazar SED sequence model. As in 
NT06, we also take into account a log-normal scatter around 
the mean relation of Lr/L 1 with a standard deviation of a p = 
0.49 in log 10 (L«/L 7 ). It should be noted that the blazar SED 
sequence is an averaged SED for groups of blazars binned by 
radio luminosity, and some scatter of Lr/L 1 is expected for 
individual blazars. 
The likelihood function is given by 



i=i 



1 dN\z u F^ u Sli) 



dzdF^dVl 



(11) 



where the subscript i is identification number of each blazar, 
A'obs the observed number of blazars, and N exp the expected 
number of the blazar detections, i.e., 



dz / dF~ / dVL 



d 3 N 



dzdFjdfl 



(12) 



The likelihood function does not depend on the normalization 
of GLF, and the normalization parameter n is determined by 
requiring A^xp = AWs- There are A^bs = 46 blazars in the sam- 
ple analyzed in NT06 and this work. 

4.2. The Best Fit Parameters 

In the first analyses, we take q as the only one free param- 
eter of the blazar GLF, with all the XLF model parameters 
fixed at the values of U03 and H05. These are hereafter called 
as U03(g) and H05(^) fits, respectively. In the second anal- 
yses, we take the faint-end slope index of XLF 71 as another 
free parameter in addition to q. These are hereafter called 
as U03(g, 71) and H05(g, 71) fits, respectively. The motiva- 
tion of these models is to investigate the effect of the uncer- 
tainty about the faint-end slope, because it may have signifi- 
cant effect on EGRB if 71 > 1 . The best-fit values for these 
four fits are shown in Table. [2] Figure [2] shows the distribu- 
tions of redshift and gamma-ray luminosity predicted by the 
best-fit models, in comparison with the EGRET data. We ob- 
tained the best-fit values of q ~ 5, meaning that the observed 
jet luminosity P (bolometric blazar luminosity) is typically 
10 5 higher than the disk X-ray luminosity. In the U03(g, 71) 
model, the characteristic AGN X-ray luminosity L\ (the break 
in XLF) corresponds to the blazar gamma-ray luminosity of 
erg s" 



We performed the Kolmogorov-Smirnov (KS) test to see 
the goodness of fits for the best-fit results of each of the 
four fits, and the chance probabilities of getting the observed 
KS deviation are shown in Table [2] Except for the H05(g) 
model, the fits are statistically acceptable. Since the best KS 
test value is obtained for the U03(g,7i) GLF model, we use 
the UQ3(q, 71) GLF model as the baseline below, when only 
one of the four fits is presented. In Fig. [3] we show the al- 
lowed regions of the model parameters in the U03(q r , 71) and 
H05(g, 71) models. The best-fit value of 71 in the U03(g, 71) 
model is in good agreement with the original value derived 
by U03. The value in the H05(<7, 71) model is slightly larger 
than the original H05 value, but the statistical significance is 
not large. The difference between the U03{q, 71) and H05(g, 
71) fits may be a result of different selections of AGNs in soft 
(H05) and hard X-ray (U03) bands, because strongly obscured 
AGNs would more easily be missed in soft X-ray bands than 
in hard X-ray bands. 

It should be noted that the results of the acceptable KS prob- 
abilities and the similar 71 values to the original XLFs give 
some support to our basic assumption of a simple proportion- 
ality between the X-ray accretion luminosity Lx and the jet 
luminosity P. As mentioned above, X-ray faint AGNs may 
be in the RIAF mode (low-Eddington ratio) rather than the 
standard disk mode, and the proportionality between Lx and 
P could be violated. Therefore our results implies that the 
EGRET data can be described well with the assumption that 
the majority of EGRET blazars are in the standard disk mode 
or high Eddington ratio phase. Fermi will probe much fainter 
blazars than EGRET, and we may see the deviation from the 
proportionality between Lx and P. 

5. THE GAMMA-RAY BACKGROUND SPECTRUM 

5.1. the EGRB Spectrum Calculation 

We calculate the EGRB spectrum by integrating our blazar 
SED sequence model in the redshift and luminosity space, us- 
ing the blazar GLF derived in $3] The spectrum of EGRB 
radiation (photon flux per unit photon energy and per stera- 
dian) is calculated as 

pL-y (/egret ,z) 



d^ie^) 
de-ydfl 



dz 



dt 



4tt J dz JL 
(1+z) L I/ [e 7 (l+z)/ft iJ ,/ , (L 7 )] 



dL 1 p 7 (L 7 ,z)(13) 



h p 

xe- T ^(z,e 7 ) , 



e 7 (l+z) 



L* = 10 48 ~" 



where e 7 is the observed gamma-ray photon energy, t the 
cosmic time, and dt/dz can be calculated by the Friedmann 
equation in the standard cosmology. We assume z max = 5 
in this calculation, but the EGRB flux is hardly dependent 
on this parameter, since the peak of GLF/XLF evolution is 
well below z ~ 5. The minimum gamma-ray luminosity 
is set at L 7 , m in = 10 43 erg s -1 as mentioned earlier. Since 
EGRET has already resolved bright gamma-ray sources, we 
should not include those sources in the EGRB calculation. 
Therefore, the maximum gamma-ray luminosity in the inte- 
gration should be L 7 (Regret j z), which is the luminosity of 
the blazar having a flux F 1 = Regret at a given z, where 
^egret = 7 x 10~ 8 photons cm" 2 s" 1 is the EGRET sensitiv- 
ity above 100 MeV. 

Very high energy photons (> 20 GeV) from high redshift 
are absorbed by the interaction with the cosmic infrared back- 
ground (CIB) radiation (Salamon & Stecker 1998; Totani & 



6 





1.8 




1.6 




1.4 




1.2 


o 




t>JQ 


1 


o 




33 


0.8 








0.6 




0.4 




0.2 








TTTTl 1 — I I I I INI 1 — I I I I INI 1 — I T 



U03(q) 

U03(q,Yi) — - 

H05(q) 

H05(q,Y x ) — ■■ 

EGRET blazars - J/': 
i'i '■ 



.1 




0.01 



0.1 
Redshift z 



1 

0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 




0.01 



0.1 

Redshift z 



III 1 1 1 1 1 III 1 1 1 1 1 III 




~ U03(q) 




- r03(q.y,) — - j 




- H05(q) fj 




" H05(q,Yi) — " // 




r EGRET blazars 

'/ ' 

'/ ,' 

if • 




~ ~ Jtf • 

A V •' 

-- J/ 








MUll ri 1— MH.LLIil*!^' 1 


1 1 f 



o 
CO 

o 



V 



0.8 
0.7 
0.6 - 
0.5 - 
0.4 
0.3 - 
0.2 - 
0.1 




1 I 1 I 1 I 1 I 

U03(q) 

U03(q,Y x ) — - 

H05(q) 
H05(q,Y x ) 
EGRET blazars 



"i — I — r 




43 



44 45 46 47 48 49 50 
logiof^yfergs" 1 ]) 



1 



0.9 - 
0.8 - 
0.7 
0.6 

0.5 H 
0.4 
0.3 
0.2 - 
0.1 - 




T" 




U03(q) 

U03(q,Y x ) — - 

H05(q) 

H05(q,Y x ) 
EGRET blazars 



43 44 45 46 47 48 49 50 
l°g 10 ( L y t er 8 s_1 ]) 



FIG. 2. — Top panels: redshift and gamma-ray luminosity (vL v at 100 MeV) distributions of EGRET blazars. The histogram is the binned EGRET data, with 
one sigma Poisson errors. The four model curves are the best-fits for the GLF models of U03(g), H05(i?), U(<?,7i), and H05(^.7i), as indicated in the figure. 
Bottom panels: the same as the top panels, but for cumulative distributions. 




FIG. 3. — Left: Solid contours show the 68%, 95%, and 99% CL likelihood contours for the LDDE model parameters [the faint-end slope index, 7j, and the 
ratio of the blazar emission power to disk X-ray luminosity, q] in the case of the U03 XLF The best-fit values (g,7i)=(4.93, 0.93) are shown by the cross. The 
dotted contours are for the parameter k, the normalization ratio of blazar GLF to AGN XLF. The k values for the contours are indicated in the figure. The shaded 
region indicates the \a error region of the original -yi value determined for AGN XLF by X-ray surveys. Right : The same as the left panel, but for the H05 XLF. 
The best-fit values (<jr,7i)=(5.35, 1.11) are shown by the cross. 
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Takeuchi 2002; Kneiske et al. 2004; Stecker, Malkan, & 
Scully 2006), and t 77 (z, e 7 ) is the optical depth of the uni- 
verse against this reaction. In this paper, we adopt the model 
of Totani & Takeuchi (2002) for CIB and r 77 . The gamma- 
ray absorption produces electron-positron pairs, and the pairs 
would scatter the cosmic microwave background (CMB) ra- 
diation to make the secondary contribution (the cascade com- 
ponent) to high energy background radiation (Aharonian et al. 
1994; Fan et al. 2004). We take into account this cascad- 
ing emission to calculate EGRB spectrum, considering only 
the first generation of created pairs, by using the same for- 
mulations in Kneiske & Mannheim (2008). In the following 
results, we will find that the amount of energy flux absorbed 
and reprocessed in intergalactic medium (IGM) is only a small 
fraction of the total EGRB energy flux, and hence the model 
dependence of CIB or the treatment of the cascading compo- 
nent does not have serious effects on our conclusions in this 
work. 

5.2. The EGRB Spectrum from Blazars 

Figure [4] shows the vF u EGRB spectrum predicted by the 
best-fit GLF model parameters obtained above. Here, we 
show the total EGRB flux as well as the contributions from 
different L 7 ranges. This EGRB spectrum is the primary spec- 
trum directly from blazars, because in this calculation the ef- 
fect of photon absorption in IGM is included but the cascade 
emission is not. The data of SMM (Watanabe et al. 1999), 
COMPTEL (Kappadath et al. 1996) and EGRET (Sreekumar 
et al. 1998; Strong et al. 2004a) are also shown. As seen in 
Figure |4] the contribution from low-luminosity blazars with 
log 1() (L 7 /erg s -1 ) = 43-45 is significant only above 1 GeV, 
because of the assumed SED sequence of blazars. The con- 
tribution from these low-luminosity blazars is larger in the 
model of H05(g,7i) than that of H05(g), because the faint- 
end slope of the GLF becomes steeper when 71 is treated as a 
free parameter. 

An important implication here is that the contribution from 
the so-called MeV blazars, which have their SED peaks at 
around MeV, is negligible in the MeV gamma-ray background 
flux, although such MeV blazars have been suspected as a 
possible origin of the MeV background. This is because MeV 
blazars have a large luminosity of log 10 (L 7 /erg s" 1 ) = 49-50 
in the blazar SED sequence, and the number density of such 
blazars is small. 

Figure [5] shows the intrinsic (the spectrum without tak- 
ing into account the absorption by CIB), absorbed (the same 
as Fig. |4), and cascading components of EGRB spectrum, 
as well as the total (absorbed+cascade) spectrum. The ab- 
sorption by CIB becomes significant at e 7 > 100 GeV, and 
the absorbed EGRB photons are converted into lower energy 
gamma-rays, with the energy flux roughly conserved. How- 
ever, the absorbed energy flux of EGRB above ~ 100 GeV is 
not significantly larger than the unabsorbed energy flux, and 
hence the increase of EGRB flux or change of the EGRB spec- 
trum at < 100 GeV by the cascading component is not a large 
effect. 

5.3. The EGRB Spectrum from Non-blazar AGNs 

In addition to blazars, we take into account non-blazar 
AGNs as the source of EGRB. Inoue et al. (2008, hereafter 
ITU08) has shown that non-blazar AGNs are a promising 
source of EGRB at ~ 1— 10 MeV. In this scenario, a nonther- 
mal power-law component extends from hard X-ray to — 10 
MeV band in AGN spectra, because of nonthermal electrons 



that is assumed to exist ubiquitously in hot coronae around 
AGN accretion disks. The existence of such nonthermal elec- 
trons is theoretically reasonable, because the hot coronae are 
the essential component in the standard picture of X-ray emis- 
sion from AGNs, and magnetic reconnection is the promising 
candidate for the heating source of the hot coronae (Laor & 
Behar 2008). Magnetic reconnections should produce non- 
thermal particles, as is well known in e.g., solar flares or the 
earth magnetosphere. Although several sources have been 
proposed as the origin of the MeV background, this model 
gives the most natural explanation for the observed MeV 
background spectrum that is a simple power-law smoothly 
connected to CXB. 

Theoretically, there is no particular reason to expect a cut- 
off of the nonthermal emission around 10 MeV, and it is well 
possible that this emission extends beyond 10 MeV with the 
same power index. Then, this component could make some 
contribution to EGRB in the EGRET energy band. The EGRB 
spectrum from the non-blazar AGNs calculated based on the 
model of ITU08 is shown in Figure [6] in comparison with the 
blazar component. The model parameters of ITU08 have been 
determined to explain the EGRB in the MeV band. There 
is some discrepancy between the MeV EGRB data of SMM 
and COMPTEL, and the reason for this is not clear. Here 
we use two model parameter sets of (r,7 tr ) = (3.5,4.4) and 
(3.8,4.4) in ITU08, where T is the power-law index of non- 
thermal electron energy distribution and 7 tl is the transition 
electron Lorentz factor above which the nonthermal compo- 
nent becomes dominant. These two parameter sets are chosen 
so that they fit to the COMPTEL and SMM data, respectively. 
The EGRB flux is dominated by blazars at the energy range 
above — 100 MeV, while the non-blazar component is domi- 
nant at the lower photon energies. 

5.4. On the Origin of EGRB 

Now we compare the EGRB data above 100 MeV with our 
model predictions. Because of the uncertainties in the EGRB 
measurements, we plotted three types of EGRB EGRET data 
in this figure (see below for the explanations of the data 
points). As can be seen in Fig. [6] the overall background 
spectrum from X-ray to 1 GeV predicted by our blazar plus 
non-blazar model is similar to the observed data. Especially, 
the EGRB prediction using the ITU08 non-blazar background 
model with T = 3.5 is in nice agreement with the observed 
data of Strong et al. (2004a) in 0.1-1 GeV. In this case, the 
predicted EGRB flux at 100 MeV can account for 80% of the 
observed flux, which is a considerably higher fraction than 
those in previous studies (Chiang & Mukherjee 1998; Miicke 
& Pohl 2000; NT06). It should be noted that the contribution 
to the EGRB flux at 100 MeV from blazars is -45 %, which 
is similar to those estimated by previous studies. The high 
fraction is due to the addition of the EGRB component from 
non-blazar AGNs. 

The prediction is still 20 % short of the observed flux, and 
the discrepancy becomes more serious when the model is 
compared with the Strong et al. (2004a) data above 1 GeV or 
with the original EGRB determination by the EGRET team 
(Sreekumar et al. 1998). Especially, the apparent excess of 
the observed EGRB flux beyond 1 GeV might be a contribu- 
tion from a completely different component, e.g., dark matter 
annihilation (e.g., Oda, Totani, & Nagashima 2005). There- 
fore we carefully discuss the possible origin of the discrep- 
ancy below. 

We should first examine the uncertainties in the model pre- 
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FIG. 4. — The blazar EGRB spectrum (energy flux per unit logarithmic photon energy and per steradian). The four panels are for the four different GLF models 
of U03(g), \]03(q, 71 ), H05(g), and HQ5(q, 71). The solid curve is for the total EGRB spectrum from all blazars, and the other curves are for a particular range 
of blazar luminosity, as indicated in the figure (L 7 in units of erg/s). The effect of absorption by CIB is included, while the reprocessed cascade component is 
not included. The observed data of SMM (Watanabe et al. 1999), COMPTEL (Kappadath et al. 1996), and EGRET [Sreekumar et al. 1998 (S98); Strong et al. 
2004a (S04a)] experiments are also shown with the symbols indicated in the figure. 
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FIG. 5. — The blazar EGRB spectrum. The two panels are for the two different GLF models of U03(g, 71 ) and H05(q. 71). Three model predictions for the 
intrinsic (no absorption by CIB), absorbed, and cascade (reprocessed emission by electrons/positrons produced in IGM) components of the EGRB spectrum are 
shown by line markings shown in the figure. The solid curve is the total flux, i.e., absorbed plus cascade components. The EGRET data are the same as those in 
Fig. [4] 
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diction. Since the GLF likelihood determination is based only 
on 46 blazars, there is a statistical uncertainty of 15 and 32% 
in the normalization of GLF and the EGRB flux at 68 and 
95% C.L., respectively. The sensitivity limit of EGRET has 
been included in the GLF likelihood analysis. The EGRET 
sensitivity to a point source changes depending on the loca- 
tion in the sky, and we used sensitivity limit inferred from 
the signal-to-noise ratios of the EGRET sources as a func- 
tion of the Galactic latitude, as in NT06. The sensitivity limit 
has a —30% scatter even after binned by the Galactic lati- 
tude. When the sensitivity limit is changed by ±30%, we 
find that the GLF normalization and EGRB flux changes by 
± —25%. Since the faint-end slopes of the (g,7i) GLF mod- 
els are 71 — 1, the faint-end cut-off luminosity L 7m i n of blazar 
GLF could also be important. We therefore repeated the cal- 
culation with L~ m i n changed by a factor of 10 from the base- 
line model, but we find that the EGRB flux changes by only 
-0.8% and -1.7% for U03(q, 71) and H05(^, 71), respec- 
tively. Therefore, our result is not sensitive to L 7 m i n . 

Next we examine the possible systematic uncertainties in 
the observational determination of EGRB. The correct mod- 
eling of the foreground, i.e., the Galactic diffuse emission is 
critical for the EGRB measurement. However, there is a well- 
known problem of "GeV anomaly", which is an excess of 
the observed diffuse flux compared with the standard theo- 
retical model of the Galactic diffuse emission (Pohl & Espos- 
ito 1998; Strong et al. 2004b; de Boer et al. 2005; Kamae, 
Abe, & Koi 2005; Stecker, Hunter, & Kniffen 2008). The dif- 
ference of the EGRB data of Strong et al. (2004a) from the 
original EGRET data (Sreekumar et al. 1998) is a result of 
modifying the model of the Galactic diffuse emission to re- 
solve the GeV anomaly, demonstrating that theoretical uncer- 
tainties in the Galactic diffuse emission could have a signifi- 
cant effect on the EGRB measurements. On the other hand, 
Stecker et al. (2008) suggested that the most likely cause of 
the GeV anomaly is a systematic error in the calibration of 
the EGRET detector at photon energies beyond 1 GeV, and 
they derived the correction factors for the flux measured by 
EGRET. If their claim is correct, the correction factors should 
be multiplied to the original EGRET measurements of EGRB 
(Sreekumar et al. 1998), and such corrected data are shown 
as S98+S08 in Fig. |6] Although the overall corrected EGRB 
flux is still higher than the model predictions, the corrected 
spectrum is similar to those of models. The marked dip at ~ 
GeV and hump at higher photon energies found in the data of 
Strong et al. (2004a) are not seen. It should also be noted that 
the EGRET EGRB data at -50 MeV is considerably higher 
than the neighboring COMPTEL data, again indicating some 
systematic uncertainties in the EGRB flux estimates. 

Based on these results and considerations, we conclude that 
most, and probably all, of the EGRB flux can be explained 
by blazars and non-blazar AGNs, with luminosity functions 
that are consistent with the EGRET blazar data and the X-ray 
AGN surveys. We must await the future Fermi data for a more 
robust conclusion about this issue. 

6. PREDICTIONS FOR THE FERMI MISSION 

6.1. Expected Number of Blazars and Non-blazar AGNs 

The left panel of Figure [7] shows the cumulative distribu- 
tion of > 100 MeV photon flux of blazars. The four GLF 
models of U03(q), U03( ? , 7l ), H05(q), and H05(tf,7i), pre- 
dict that about 640, 930, 470, and 1200 blazars should be de- 
tected by Fermi, respectively, where we assumed the Fermi 
sensitivity to be F\i m = 2 x 10~ 9 photons cm" 2 s" 1 above 100 



MeV. This is a standard sensitivity often used in the literature 
(e.g., Oh 2001; NT06; Venters & Pavlidou 2007), and close to 

— 5<T limit for a point source at high Galactic latitude after a 1 
year survey. See also the official Science Requirement Docu- 
ment of Fermi available at http://fermi.gsfc.nasa.gov/science/ 
for the Fermi performance. Here, we calculated simply all the 
blazars without taking into account the probability of identi- 
fication by follow-up observations in other wavebands. This 
is why the number of bright blazars is larger than the EGRET 
observation in Fig. |7j we have taken into account the proba- 
bility of blazar identification by detecting a radio counterpart 
in the likelihood analysis of the EGRET blazars. The model of 
Stecker & Salamon (1996), and the PLE and LDDE models 
in Narumoto & Totani (2006) predicted -10000, 5400, and 
3000 blazars for the same sensitivity limit, respectively. It is 
remarkable that the GLF models in this work predict signifi- 
cantly smaller numbers of blazars than previous studies. This 
is probably because of the SED model newly used here; the IC 
component has its broad SED peak at around 100 MeV, and 
the flux at lower and higher energy band is relatively lower 
than the case of power-law SED, because of the curvature of 
the spectrum. 

The right panel of Figure|7]shows the differential flux distri- 
bution of gamma-ray blazars multiplied by flux, showing the 
contribution to the EGRB per unit logarithmic flux interval. 
From this plot one can estimate how much fraction of EGRB 
will be resolved by the Fermi mission. The peak of the contri- 
bution to EGRB occurs at a flux much brighter than the Fermi 
limit, meaning that EGRB from blazars should practically be 
resolved into discrete sources. We find that the fraction of 
EGRB flux that should be resolved is 99, 98, 100, and 100% 
against the total blazar EGRB flux in the four GLF models 
of U03(^), U03(4, 7 i), H05(^), and H05(g,7i), respectively. 
However, as seen in Fig. [6] there is a considerable contribu- 
tion from non-blazar AGNs to the EGRB flux at 100 MeV, 
and our next interest is how much fraction of EGRB from 
non-blazars will be resolved by Fermi. 

The left panel of Figure [8] shows the cumulative flux distri- 
bution of blazars in the U03(q,~f\) model, non-blazar AGNs, 
and the total of the two. The expected number of non-blazar 
AGNs detectable by EGRET is much smaller than unity in 
all sky, and it is consistent with the fact that almost all extra- 
galactic gamma-ray sources detected by EGRET are blazars. 
However, about 1-30 non-blazar AGNs would be detected 
by Fermi, by their soft nonthermal emission from nonther- 
mal electrons in coronae of accretion disks, giving an in- 
teresting test for our model. In fact, we can estimate the 
> 100 MeV gamma-ray flux from the observed hard X-ray 
flux of NGC 4151 (the brightest Seyfert galaxy in all sky, 
Sazonov et al. 2007) using the ITU08 model, and it becomes 

- 3.3 x 10" 8 and 4.1 x 10" 9 photons cm" 2 s" 1 for T = 3.5 
and 3.8, respectively. Therefore NGC 4151 is marginally de- 
tectable when r = 3.8, and easily detectable with T = 3.5, 
which are nicely consistent with the predicted source counts 
of non-blazar AGNs detectable by Fermi. 

On the other hand, the right panel of Figure ^indicates that 
EGRB from non-blazar AGNs is hardly resolved into discrete 
sources by Fermi. The predicted fraction of non-blazar EGRB 
flux resolved by Fermi is 0.03 and 0.004 % for the models 
with T = 3.5 and 3.8, respectively, against the total non-blazar 
EGRB flux. These results are in sharp contrast to those for 
blazars, although the contributions to the total EGRB by the 
two populations are comparable at around 100 MeV. This is 
because of a large difference between typical absolute lumi- 
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FIG. 6. — The EGRB spectrum from non-blazar AGNs and blazars. The two panels are for the two different blazar GLF models of U03(g,7i) (left) and 
H05(ir/,7i) (right). The model curves of the blazar component (absorbed+cascade), non-blazar AGN component, and the total of the two populations are shown. 
Note that two models are plotted for the non-blazar component with different values of T (see the line-markings indicated in the figure). The observed data of 
HEAO-1 (Gruber et al. 1999) and Swift/BAT (Ajello et al. 2008) are shown. We also plot a new EGRET data denoted as "S98+S08", which is the original 
EGRET determination of Sreekumar et al. (1998) corrected by the correction factors proposed by Stecker et al. (2008, S08). The other data are the same as those 
in Fig. |4] 
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FIG. 7. — Left: the cumulative flux distribution of blazars. The four model curves correspond to the four different GLF models of U03(g), H05(g), U(<j,7i), 
and H05(g. 71 ) are shown. The thin solid line shows the observed distribution of EGRET blazars. The detection limits of EGRET and Fermi are also shown. 
Right: the same as the left panel, but showing differential flux distribution multiplied by flux F-,, to show the contribution to EGRB per logarithmic flux interval. 
The solid squares are the EGRET data with Poisson errors. 



nosities of blazars and non-blazar AGNs. The typical lumi- 
nosity at 100 MeV of blazars and non-blazars is vL v = 10 48 
and 10 42 5 erg/s, which are corresponding to L\ in the AGN 
XLF. If the evolution is similar for the two population, the 
characteristic redshift or distance to the sources having main 
contribution to EGRB should be similar. Therefore, a source 
population with smaller characteristic luminosity is more dif- 
ficult to resolve into discrete sources when the flux sensitivity 
is fixed. 

The above results mean that a considerable part of EGRB at 
around 100 MeV will remain unresolved even with the Fermi 
sensitivity, because there is a considerable contribution from 
non-blazar AGNs to EGRB at ~ 100 MeV. However, the con- 
tribution from non-blazars should rapidly decrease with in- 
creasing photon energy, and almost all of the total EGRB 
flux at > 1 GeV should be resolved into discrete blazars by 
Fermi, if there is no significant source contributing to EGRB 



other than blazars. It should be noted that these predictions 
can be tested by Fermi relatively easily, without follow-up or 
cross check at other wavebands, once measurements of source 
counts and EGRB flux have been done by the Fermi data. 
Therefore this gives a simple and clear test for the theoreti- 
cal model presented here. 

6.2. Redshift and Luminosity Distribution 

Figure [9] shows the redshift and luminosity distributions of 
blazars that will be detected by Fermi, where we have again 
set the Fermi sensitivity limit as Fn m = 2 x 10~ 9 photons cm" 2 
s" 1 above 100 MeV. Since we normalized the total number of 
detectable blazars, only the shapes of distribution should be 
compared. 

Fermi will detect much fainter blazars than EGRET did. 
Since 71 > 1, the H05(g, 71) model predicts more faint blazars 
than other models. The redshift distribution of EGRET 
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FIG. 8. — The same as Fig. [7] but showing non-blazar AGNs as well, in addition to blazars. The two models of non-blazar AGNs with different values of T are 
shown. The total of blazar and non-blazar counts is also shown. 
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blazars has already extended to cosmologically significant 
range of z > 1 , and the normalized distribution does not signif- 
icantly change by the Fermi sensitivity. However, the absolute 
number of high-z blazars will be increased, because the total 
number of blazars will drastically be increased by Fermi from 
that of EGRET blazars. 

There are some differences in the predictions by the four 
different GLF models. Since the Fermi will dramatically in- 
crease the statistics by a large number of detected blazars, a 
detailed comparison between the GLF models and the Fermi 
data will provide us a quantitative view of the evolutionary 
nature of blazar GLF, once the redshifts of Fermi blazars are 
measured. The comparison between the blazar GLF and other 
AGN LFs in different wavebands (e.g., X-ray, optical, and ra- 
dio) will give us important information about the jet activity 
evolution of AGNs, in comparison with the mass accretion 
history onto SMBH that can be probed by disk luminosity 
function. 

7. DISCUSSION 
7.1. Comparison with previous studies 



There are four important points in the modeling of the 
blazar GLF. They are (1) the treatment of blazar population 
(e.g., BL Lac/FSRQ separation), (2) modeling of blazar SED, 
(3) the functional form of the GLF and its cosmological evo- 
lution, and (4) observational inputs to determine the GLF pa- 
rameters. Here we compare our model with those of previ- 
ous studies about these points. Our results about the expected 
number of Fermi blazars (—1000 blazars to the sensitivity 
limit of film = 2 x 10" 9 photons cm" 2 s" 1 above 100 MeV) 
and EGRB (~ 45% contribution by blazars at 100 MeV) will 
also be compared with those of the previous studies. 

Stecker & Salamon (1996; hereafter SS96) treated the 
blazar as a single population, with a single power-law SED. 
They used the blazar radio luminosity function of Dunlop & 
Peacock (1990) with pure luminosity evolution (PLE) as the 
blazar GLF model. The gamma-ray/radio flux ratio was esti- 
mated by those of typical blazars, and their blazar GLF was 
not statistically compared with the EGRET blazar data of flux 
and redshift distribution. They found that blazars can explain 
100% of EGRB, and -10000 blazars are expected to be de- 
tected by Fermi. 
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FIG. 10. — The blazar SED models in our model in comparison with those in previous studies. BL Lacs and FSRQs are separately plotted in the left and right 
panels, respectively. SEDs of our model are those having L 7 = 10 44 ' 5 and 10 45 6 ergs _1 for BL Lac and FSRQ, respectively. The MP00 curves are the EIC+SSC 
instantaneous injection model from their Fig 1. NT06 used a single power-law with the same index both for BL Lacs and FSRQs. The G06 curve (the same for 
BL Lac and FSRQ) is the SED with the synchrotron peak of ^ syn = 10 13 8 Hz used in their study. D07 used single power-law SEDs but with different indices 
for BL Lacs and FSRQs. The KM08 curve is the case 3 SED used as a model of HBL SED. The normalizations of these SEDs are from the original models for 
MP00 and D07, and otherwise they are fixed at L 7 = 10 44 5 and 10 46 6 ergs~' for BL Lac and FSRQ, respectively, where L 7 is vL v at 100 MeV. 



Chiang & Mukherjee (1998; hereafter CM98) treated the 
blazar as a single population, with a single power-law SED. 
CM98 assumed a single power-law with a cut-off luminosity 
^rnin for the luminosity distribution of GLF, and adopted PLE 
for evolution. In such a modeling, the faint-end slope is dif- 
ficult to estimate accurately, and hence the uncertainty about 
EGRB flux becomes large. They constrained the blazar GLF 
parameters from a statistical comparison with the EGRET 
blazar data. They found that blazars can explain only 25% 
of EGRB. 

Mucke & Pohl (2000; hereafter MP00) treated BL Lacs and 
FSRQs separately. They used a theoretical model of the blazar 
SED of Dermer & Schlickeiser (1993), and their GLF model 
is described in terms of the distribution of electron energy in- 
jection. The energy injection distribution is modeled accord- 
ing to the radio luminosity function of Fanaroff-Riley (FR) I 
and II galaxies of Urry et al. (1991) and Padovani & Urry 
(1992) for BL Lacs and FSRQs, respectively. These LFs as- 
sume PLE. The GLF model parameters are then determined 
by the flux and redshift distributions of EGRET blazars. They 
found that blazars can explain 20-80% of the EGRB depend- 
ing on the model parameters. 

Narumoto & Totani (2006; hereafter NT06) treated the 
blazar as a single population, with a single power-law SED. 
They examined both the PLE and LDDE model for GLE evo- 
lution, which are based on the radio LF of blazars and XLF 
of AGNs, respectively. The GLF model parameters are de- 
termined by the fits to the flux and redshift distributions of 



EGRET blazars. They found that LDDE model fits to EGRET 
data better than the PLE model. They found that blazars can 
explain 25-50% of the EGRB, and ^3000 blazars are ex- 
pected to be detected by Fermi. 

Giommi et al. (2006; hereafter G06) estimated EGRB 
flux and spectrum based on the radio luminosity function of 
blazars. They first calculated the cosmic radio background in- 
tensity, and then converted it to the EGRB spectrum by using a 
typical SSC SED model. Therefore, their model is not a GLF 
model for individual blazars, and not compared statistically 
with the EGRET data of blazar flux and redshift distribution. 
Their model overproduced the observed EGRB spectrum. 

Dermer (2007; hereafter D07) treated BL Lac and FSRQ 
separately, using single power-law SEDs for the two popu- 
lations with different spectral indices. The variation of the 
absolute luminosity was not considered within these popula- 
tions, with the single luminosity of L 7 = 2.5 x 10 44 erg s" 1 
and 10 46 erg s" 1 for BL Lacs and FSRQs, respectively. There- 
fore, this is not a model of GLF, but D07 rather considered 
various models of number density evolution, and the model 
parameters of the number density evolution are constrained 
by the redshift and flux distribution of EGRET blazars. They 
found that blazars can accountfor 12-19% of EGRB at 1 GeV, 
and their model predicts ~ 1000 blazars detectable by Fermi, 
which is close to our result by a very different approach. 

Kneiske & Mannheim (2008; hereafter KM08) considered 
the contribution to EBL from high energy cutoff BL Lacs 
(HBLs). The spectrum was assumed to be a broken power-law 
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with the peak energy around 1 TeV for HBLs. They employed 
the XLF of HBL which shows no evolution (Beckmann et al. 
2003), and a statistical comparison with the EGRET blazar 
data was not made. They found that blazars can explain 
—80% of EGRB by including cascading emission (see also 

E2>. 

For comparison, we show the SEDs of our model cor- 
responding to the typical luminosities of BL Lacs (L 7 ~ 
10 44 - 5 erg s" 1 ) and FSRQs (L 7 - 10 46 6 erg s" 1 ), as well as 
those used in previous studies, in Fig. [10] 

It should be noted that, in the studies that constrained the 
blazar GLF parameters by a statistical comparison with the 
EGRET data, only CM98, NT06, and this work take into ac- 
count the probability of identification of gamma-ray blazars 
by finding radio counterparts. CM98 assumed that there is 
no correlation between gamma-ray and radio luminosities of 
blazars. The assumption of no correlation at all over a wide 
range of gamma-ray and radio luminosities induces some in- 
consistencies (see discussion in Stecker & Salamon 2001), 
and it is physically reasonable to expect some level of cor- 
relation between the two from the viewpoint of the standard 
synchrotron and inverse Compton emission model of blazars. 
In NT06 and this work, we introduced this correlation with 
a log-normal scatter in the gamma/radio flux ratio, which is 
consistent with the observation. 

7.2. Cascade of Very High Energy Gamma-Rays and EGRB 

KM08 estimated that HBL would account for —20% of the 
EGRB by cascading emission produced by absorption of very 
high energy gamma-rays in IGM. However, in our estima- 
tion the cascade emission would not significantly contribute 
to EGRB (see Fig[5]l, and the reason for this discrepancy is as 
follows. In KM08, the minimum and maximum luminosities 
of HBL are set as L m j n = 10 43 erg s" 1 and L max = 10 47 erg s -1 , 
respectively (in vL v at 0.3 TeV). Their SED model is a 
broken-power law, which is independent of luminosity and 
obtained from fitting to the TeV observations of Mkn421, 
PKS2005-489, and PKS2155-304. However the luminosities 
of these three at 0.3 TeV are about 10 44 erg s" 1 , which is rela- 
tively low in the range of HBL luminosity assumed by KM08. 
The use of the same SED up to L max = 10 47 erg s" 1 may over- 
predict the flux in very high energy band, compared with the 
case of the SED sequence. Therefore we conclude that the 
contribution to EGRB from the cascading component is not 
significant provided that the SED sequence is valid. 

7.3. Effects of Flaring Activities 

It is well known that blazars show flaring activities. For 
example, TeV Cherenkov telescopes H.E.S.S. and MAGIC 
have recently observed minute-scale variability in TeV ranges 
(Aharonian et al. 2007; Albert et al. 2007). The origin of such 
variability has been discussed in many papers (Begelman et al. 
2008; Ghisellini & Tavecchio 2008a; Ghisellini et al. 2008; 
Katarzyhski et al. 2008; Mastichiadis & Moraitis 2008; Gian- 
nios et al. 2009), but it is still rather uncertain. Therefore it is 
difficult to model the flaring activity quantitatively and incor- 
porate in our analysis. We used the time-averaged gamma-ray 
flux of blazars in the third EGRET catalog, which is the mean 
of several viewing periods during the total EGRET observa- 
tion period (see §2.1 of NT06 for more details). Typical dura- 
tion of one viewing period is 2-3 weeks, and flaring activities 
within this scale have already been taken into account. Flar- 
ing activities that might have happened at epochs out of the 



viewing periods are not taken into account, but still the flaring 
effect should have been incorporated to some extent by the use 
of time averaged flux. More observational information about 
blazar flares (light curves and frequency, etc.) is required to 
examine this issue in more detail, and future Fermi data will 
be useful. 

8. CONCLUSIONS 

In this paper, we constructed a new model of the blazar 
gamma-ray luminosity function (GLF), taking into account 
the blazar SED sequence and the LDDE luminosity func- 
tion inferred from X-ray observations of AGNs. An im- 
plicit assumption here is that the jet activity of AGNs is as- 
sociated with the accretion activity, and hence the blazar lu- 
minosity function has a similar evolution to that of X-ray 
AGNs. The GLF model parameters are constrained by care- 
fully fitting to the observed flux and redshift distribution of the 
EGRET blazars. By this model, for the first time, we can pre- 
dict the spectrum of the extragalactic gamma-ray background 
(EGRB) in a non-trivial way, rather than assuming a simple 
functional form such as power-law spectra. 

The absorption of gamma-rays in IGM by interaction with 
the cosmic infrared background (CIB), and the reprocessed 
cascade emission from electrons/positrons produced in IGM 
are also taken into account in the EGRB calculation. We 
found that, provided that the blazar SED sequence is valid, 
the amount of EGRB energy flux absorbed by CIB interaction 
(at photon energy > 10 GeV) is not large and the EGRB spec- 
trum at lower photon energy bands is not significantly affected 
by the secondary cascade emission. 

The contribution from non-blazar AGNs to EGRB is also 
considered, by using the nonthermal coronal electron model 
of Inoue et al. (2008). This model is a natural extension of the 
standard X-ray spectral model of AGN emission from the disk 
and corona region, and predicts the nonthermal emission ex- 
tending from MeV to higher energy band with a steep power- 
law. This model gives a natural explanation for the observed 
cosmic MeV background, and we examined whether X-ray to 
GeV gamma-ray background radiation can be accounted for 
by the two population model including blazars and non-blazar 
AGNs. 

Our model predicts that the EGRB flux from blazars at 100 
MeV is only about 45 % of the observed EGRB flux of Strong 
et al. (2004a). However, it is possible to account for more 
than 80% of the observed EGRB flux by the sum of blazars 
and non-blazar AGNs. The predicted spectrum is also in good 
agreement with the observed data from X-ray to — 1 GeV, 
indicating that the EGRB below 1 GeV can mostly be ex- 
plained by our two population model. The two components 
have a comparable contribution to EGRB at —100 MeV, and 
the blazar component is dominant at the higher energy range 
while the non-blazar component at the lower. Therefore, the 
higher EGRB flux at 100 MeV found by this work than many 
of previous studies is not by change in the EGRB flux from 
blazars, but by adding non-blazar AGNs. 

The EGRB spectrum of Strong et al. (2004a) shows a hump 
beyond 3 GeV, which cannot be explained by our model. 
The EGRB flux determined originally by the EGRET team 
(Sreekumar et al. 1998) is higher than our model, and the ob- 
served spectrum is harder. However, measurements of EGRB 
suffer from the uncertainties in the theoretical modeling of the 
foreground diffuse Galactic emission, and perhaps from sys- 
tematic errors in the calibration of the EGRET detector above 
1 GeV (Stecker et al. 2008). We must await the Fermi data to 
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examine these issues and whether an additional, completely 
different source is necessary to explain all of EGRB. We con- 
clude that most (at least more than 50%), and probably all, of 
the observed high-energy background radiation from X-ray to 
1 GeV spanning 6 orders of magnitude in photon energy can 
be accounted for by AGNs including blazars. 

We predicted the flux distribution of blazars by our model, 
and we found that 600-1200 blazars in all sky should be de- 
tected by Fermi, assuming a sensitivity limit of F\\ m = 2 x 10~ 9 
photons cm" 2 s" 1 above 100 MeV. This number is signifi- 
cantly lower than the predictions by many of the previous 
studies. The expected number of the non-blazar AGNs at 
100 MeV band is much smaller than unity in all sky at the 
EGRET sensitivity, but about 1-30 of such population are ex- 
pected be detected by Fermi. Fermi should resolve almost 
all of the EGRB flux from blazars into discrete blazars, with 
percentages of >99%. On the other hand, less than 0.1% 
of the EGRB flux from non-blazar AGNs can be resolved 
into discrete sources. Therefore, we have a clear prediction: 
Fermi will resolve almost all of the EGRB flux into discrete 
sources at photon energies > 1 GeV where blazars are dom- 
inant, while a significant fraction of the EGRB flux will re- 
main unresolved in the low energy band of < 100 MeV where 
non-blazar AGNs have a significant contribution. This predic- 
tion can easily be tested, only with the source counts and the 
EGRB estimates by Fermi data. 



We also predicted the redshift and absolute luminosity dis- 
tributions for the Fermi blazars. Future Fermi data set with 
measured redshifts will enable us to get a quantitative mea- 
surement of the blazar GLF and its evolution with a much 
larger statistics than EGRET. A direct comparison between 
blazar GLF and X-ray AGN LF could be done, making it pos- 
sible to discuss the relation between the cosmic histories of 
jet activity and accretion activity of AGNs. Our model will 
be able to serve as a guide in such studies. For example, if 
the jet luminosity is proportional to the mass accretion rate, 
we expect a larger jet/disk ratio than assumed in this work 
for AGNs with low Eddington ratio, because X-ray luminos- 
ity scales as oc m 2 in the RIAF regime. Then, we may expect 
a different behavior of blazar GLF from that of AGN XLF at 
faint luminosity range. Such an evidence is already seen from 
our fit to the EGRET data (see Fig. and it will be tested 
more quantitatively by the Fermi data soon. 
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APPENDIX 

THE BLAZAR SED SEQUENCE FORMULATIONS 

We define ip(x) = log 10 [^L„/(erg s" 1 )] with x = log 10 (V/Hz) (v in rest-frame). The empirical SED sequence model of blazars 
is the sum of the synchrotron [ip s (x)] and IC [ip c (x)] emissions, i.e., 

ip(x) = log 10 [10^ w + 10*« w ] . (Al) 

Each of the two components is described by a combination of a linear and a parabolic functions at low and high photon frequen- 
cies, and we define x tr s and x tr , c as the linear-parabolic transition frequencies for the synchrotron and IC components, respectively. 
We take tp R = log 10 [L R /(erg s" 1 )] as a reference of a blazar luminosity, where L R is vL v luminosity in the radio band (v R =5 GHz 
or x R = 9.698). The linear part of the synchrotron component is described as follows: 

ijj s (x) = ip s i(x) = (l-a s )(x-x R ) + 4> R (x < x a . s ) , (A2) 

where a s = 0.2 is the energy flux index (i.e., L v oc v~ a '). The parabolic part of the synchrotron component is characterized by the 
vL v peak frequency v s (or corresponding x s ), as 

ip s (x) = tp s2 (x) = -[(x-x s )/a] 2 + tp s .p (x > x a . s ) , (A3) 

where a is a parameter that controls the width of the parabolic function. The peak luminosity tp S:P of the synchrotron component is 
determined if x tr , s , x s , and a are given, by requiring the continuity between the linear and parabolic parts, i.e., ip s i(xtr, s ) = ip S 2(xti,s)- 
The result is: 

ip S)P = (l-a s )(x l[ , s -x R ) + ip R + ^ Xas ^ . (A4) 

The linear part of the IC component (roughly in the hard X-ray band) is defined as: 

ij) c (x) = ip c i(x) = (l-a c )(x-x x ) + ipx (x < xtr, c ) , (A5) 

where the power index is a c = 0.6 (different from the synchrotron component), and tp x is the luminosity of IC component at the 
reference frequency vx = 2.42 x 10 17 Hz = 1 keV/h p , orxx = 17.383. The parabolic part of the IC component is characterized by 
the vL v peak frequency v c (or corresponding to x c ) with the same width parameter a as the synchrotron component, i.e., 

1p c (x) = 1p c2 {x) = -[(x-X c )/a] 2 + lp C! p (X > Xtr, c ) , (A6) 

where ip CjP is the peak luminosity of the IC component. The linear-parabolic transition frequency x^ is determined by requiring 
continuity, ip c \(x a .c) = ipdixa.c) and the result is: 



where 



x tr,c — ^ ) (A7) 

C = cr 2 (l-a c )-2x c , (A8) 

r] = x 2 . + <7 2 [if)x-x x (l-a c )-ip c . p ] . (A9) 

Then, the SED sequence is determined whenx^, x„ cr, ip x , x c , and ip c , p are given as functions of ip R . First we setx tr s =log 10 (5 x 
10 10 ) and lO^ - *' = v c jv s = 5 x 10 8 , for all luminosity range of blazars. To get a good fit in all luminosity range, we divide the 
luminosity into two ranges of ip R < 43 and ip R > 43. The parameters relevant for the synchrotron component in the ip R < 43 
range are determined as follows. The synchrotron peak frequency is determined by 

x s = -ai>R -43) + 14.47, (A10) 

with £ = 0.88. The width parameter a is determined by the following scaling law: 

a = 0.089 lx s + 1.78 . (All) 

In this case, the connection atx = x tr s is not smooth (i.e., not differentiable), but the discontinuity of the derivative is not significant. 
In the range of ip R > 43, we change the synchrotron peak frequency x s by adopting £ = 0.4, and the width parameter a is 
determined so that the connection at x tr s is differentiable, i.e., 

a = [2(x s -x tr , s )/(l-« s )] 1 / 2 . (A12) 
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For the IC component parameters, the IC peak luminosity ip CtP is assumed to be the same as the synchrotron peak, i.e., ip CtP = 
ip s ,p, at tpR < 43, while it is determined by the following formula at ipR > 43: 



where (3 = 1.77, 6 = 0.718, and e = 45.3. Still, the normalization of the IC linear part, ipx, remains to be determined. To get a 
good fit to the X-ray data, we define the dependence of this parameter on ip R as follows: 



with ipx,43 = 43.17. The parameter ip x is kept constant at ip R > 46.68, because ip c i(x) = VtaW does not have solution and hence 
Xtrx cannot be obtained, when we extrapolate the ip x formula at ip R < 46.68 into ip R > 46.68. At tp R = 46.68, il> CtP becomes 49.81, 
which is brighter than the maximum L 7 of the EGRET blazars, and is larger than the characteristic L* corresponding to the break 
X-ray luminosity L* x in the AGN XLF by a factor of ~ 60. Therefore, the treatment of ipx at such large luminosity range is not 
important when we consider the overall properties of blazars and EGRB. The linear-parabolic connection of the IC component at 
Xtr c is not as smooth as that of the synchrotron component at x^s, but this is inevitable to fit the SED in the X-ray region. 

It should be noted that the formulations of x s , a, ip C:P , and ipx are continuous at ip R = 43 and 46.68. Therefore the SED in this 
formula always changes continuously with tp R in the whole luminosity range, which was the motivation of constructing this new 
formulation. 



(A13) 




Wr < 43) 

(43 < ip R < 46.68) 

(46.68 < ip R ) 



(A14) 



